Reference sites

Packages

pacman::p_load(
  dplyr,
  tidyverse,
  readxl,
  vegan)
# setwd("D:/PhD/01_Data/01_Baseline/Model_1") # Windows
# setwd("/Volumes/Yiming_Wang/PhD/01_Data/01_Baseline/Model_1") # Mac
# setwd("/Users/godric/Desktop/Fut2 review comments/Revision") #Yiming's macbookpro

setwd("/Volumes/Yiming_Wang/PhD/15_Manuscript/Submission/07_ISME/Review comments/Revision data/26 April 2024")#Yiming's imac pro
df_all <- read_excel("Baseline_S33_taxa_formatted.xlsx")

# Subgroup by variables
#30436d WT Male
#9ea5c2 WT Female
#654e3a KO Male
#bcaba6 KO Female

#3C5488B2 WT
#7E6148B2 KO
#D98594FF Female
#94C5CCFF Male

Alpha diversity_visualization and analysis

Plot-all mice

# load package
pacman::p_load(
  ggplot2,ggsci, ggthemes,
)

set.seed(123)

# Level orders
level_order_genotype <- c("WT","KO")
level_order_sex <- c("Male","Female")
level_order_genotypesex <- c("WT Male","WT Female","KO Male","KO Female")

# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal
## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"
library("scales")
show_col(mypal) #  #3C5488B2(WT) and  #E64B35B2 (KO)

# Genotype
## FDP Diversity
Genotype_FDP <- ggplot (data=df_all, mapping = aes(y=Diversity_FPD, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(4,24),breaks = seq(4,24,5))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=24,
           label="ns",hjust=0.5, size=5.5)


Genotype_FDP_V <- Genotype_FDP + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)

## Shannon Diversity
Genotype_shannon <- ggplot (data=df_all, mapping = aes(y=Diversity_Shannon, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(5,8),breaks = seq(5,8,1))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=8,
           label="ns",hjust=0.5, size=5.5)


Genotype_shannon_V <- Genotype_shannon + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)

## Richness (Chao)

Genotype_chao <- ggplot (data=df_all, mapping = aes(y=Richness_Chaos, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Chao richness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(50,250),breaks = seq(50,250,50))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=250,
           label="ns",hjust=0.5, size=5.5)


Genotype_chao_V <- Genotype_chao + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)


## Pielou evenness
Genotype_pielou <- ggplot (data=df_all, mapping = aes(y=Evenness_Pielou, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(0.8,1.0),breaks = seq(0.8,1.0,0.05))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=1.0,
           label="ns",hjust=0.5, size=5.5)


Genotype_pielou_V <- Genotype_pielou + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)






# Print
Genotype_FDP

Genotype_shannon

Genotype_chao

Genotype_pielou

Plot-male mice only

# Level orders
level_order_genotype <- c("WT","KO")
level_order_sex <- c("Male","Female")
level_order_genotypesex <- c("WT Male","WT Female","KO Male","KO Female")

# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal
## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"
library("scales")
show_col(mypal) #  #3C5488B2(WT) and  #E64B35B2 (KO)

df_male <- df_all %>% 
  filter(Sex == "Male")

## FDP Diversity
Genotype_FDP_male <- ggplot (data=df_male, mapping = aes(y=Diversity_FPD, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(4,24),breaks = seq(4,24,5))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=24,
           label="ns",hjust=0.5, size=5.5)


Genotype_FDP_V <- Genotype_FDP + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)

## Shannon Diversity
Genotype_shannon_male <- ggplot (data=df_male, mapping = aes(y=Diversity_Shannon, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(5,8),breaks = seq(5,8,1))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=8,
           label="ns",hjust=0.5, size=5.5)


Genotype_shannon_V <- Genotype_shannon + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)

## Richness (Chao)
Genotype_chao_male <- ggplot (data=df_male, mapping = aes(y=Richness_Chaos, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Chao richness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(50,250),breaks = seq(50,250,50))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=250,
           label="ns",hjust=0.5, size=5.5)


Genotype_chao_V <- Genotype_chao + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)


## Pielou evenness
Genotype_pielou_male <- ggplot (data=df_male, mapping = aes(y=Evenness_Pielou, x=factor(Genotype, level=level_order_genotype))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(0.8,1.0),breaks = seq(0.8,1.0,0.05))+
  scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) + 
  annotate("text",x=1.5,y=1.0,
           label="ns",hjust=0.5, size=5.5)


Genotype_pielou_V <- Genotype_pielou + geom_violin(alpha=0, width=0.9,
                                             position=position_dodge(width = 0.8),
                                             size=0.75)






# Print
Genotype_FDP_male

Genotype_shannon_male

Genotype_chao_male

Genotype_pielou_male

Beta diversity_results and visualisation

Datashape

# places that are not satisfied
  ## 1) Permanova results table
  ## 2) Figure A and B did not align because of legend (Change legend size, position to bottom, change P value position)


#library packages
pacman::p_load(
  data.table,
  vegan,
  knitr,
  EcolUtils
)

# Shape data
# Level orders
df_all$Genotype <- factor (df_all$Genotype, levels = c("WT", "KO"))
df_all$Sex <- factor(df_all$Sex, levels = c("Male", "Female"))
df_all$GenotypeSex <- factor(df_all$GenotypeSex, levels = c("WT Male", "WT Female", "KO Male","KO Female"))


# subset dateframe
df_abundance <- subset(df_all, select = -c(1:17))
df_metadata <- subset(df_all, select = c(1:17))

Bray-curtis

set.seed(123)
 # Check the stress and choose the optimal number of dimensions
dist_bray <- vegdist(df_abundance,  method = "bray")

# NMDS.scree(dist)
  ## For a good representation of your data, the acceptable stress value: <0.2
  ## Considering our stress value is 0.13, which is good enough for representation of microbiome composition, so k=2 is chosen for NMDS plot
  ## The stress value does not affect PERMANOVA results, only tells you whether 2D (k=2) is good for composition demonstration

# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)
## Run 0 stress 0.1413421 
## Run 1 stress 0.1415037 
## ... Procrustes: rmse 0.004186663  max resid 0.02572238 
## Run 2 stress 0.1416862 
## ... Procrustes: rmse 0.007319768  max resid 0.04434329 
## Run 3 stress 0.1413421 
## ... Procrustes: rmse 2.594225e-05  max resid 0.0001363623 
## ... Similar to previous best
## Run 4 stress 0.1413421 
## ... Procrustes: rmse 5.953767e-06  max resid 3.293415e-05 
## ... Similar to previous best
## Run 5 stress 0.1413421 
## ... Procrustes: rmse 3.72446e-05  max resid 0.0002456958 
## ... Similar to previous best
## Run 6 stress 0.14154 
## ... Procrustes: rmse 0.005891884  max resid 0.04412429 
## Run 7 stress 0.1415037 
## ... Procrustes: rmse 0.004208291  max resid 0.02583719 
## Run 8 stress 0.1416861 
## ... Procrustes: rmse 0.007331071  max resid 0.04434624 
## Run 9 stress 0.1844246 
## Run 10 stress 0.1615297 
## Run 11 stress 0.161271 
## Run 12 stress 0.1413422 
## ... Procrustes: rmse 8.19954e-05  max resid 0.0004571962 
## ... Similar to previous best
## Run 13 stress 0.1413421 
## ... New best solution
## ... Procrustes: rmse 4.073193e-06  max resid 2.546618e-05 
## ... Similar to previous best
## Run 14 stress 0.1612711 
## Run 15 stress 0.1413421 
## ... Procrustes: rmse 1.713322e-05  max resid 0.0001087092 
## ... Similar to previous best
## Run 16 stress 0.1413421 
## ... Procrustes: rmse 2.929669e-05  max resid 0.0001228787 
## ... Similar to previous best
## Run 17 stress 0.1416861 
## ... Procrustes: rmse 0.007329494  max resid 0.04434847 
## Run 18 stress 0.1413421 
## ... Procrustes: rmse 1.40533e-05  max resid 9.236742e-05 
## ... Similar to previous best
## Run 19 stress 0.1413423 
## ... Procrustes: rmse 3.979513e-05  max resid 0.000235021 
## ... Similar to previous best
## Run 20 stress 0.1415037 
## ... Procrustes: rmse 0.004176246  max resid 0.02566424 
## *** Solution reached
NMDS$stress
## [1] 0.1413421
# PERMANOVA and Pair-wise PERMANOVA

set.seed(123)
## Genotype and sex are crossed variable, cage are nested in genotype and sex: https://stats.stackexchange.com/questions/350462/can-you-perform-a-permanova-analysis-on-nested-data
adonis2(df_abundance ~ (Genotype+Sex)/Cage, data=df_metadata, permutations=9999, method="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ (Genotype + Sex)/Cage, data = df_metadata, permutations = 9999, method = "bray")
##                   Df SumOfSqs      R2      F Pr(>F)    
## Genotype           1   0.1057 0.02780 2.7448 0.0280 *  
## Sex                1   0.0824 0.02166 2.1391 0.0721 .  
## Genotype:Sex:Cage 16   1.8816 0.49481 3.0537 0.0001 ***
## Residual          45   1.7330 0.45573                  
## Total             63   3.8027 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Genotype x Sex effect
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray") 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
##             Df SumOfSqs      R2      F Pr(>F)   
## GenotypeSex  3   0.4048 0.10645 2.3827 0.0068 **
## Residual    60   3.3979 0.89355                 
## Total       63   3.8027 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df_permanova_pairwise <- adonis.pair(vegdist(df_all[,18:128],  method = "bray"), df_all$GenotypeSex, nper = 9999, corr.method = "fdr")
kable(df_permanova_pairwise)
combination SumsOfSqs MeanSqs F.Model R2 P.value P.value.corrected
WT Male <-> WT Female 0.0700924 0.0700924 1.2445858 0.0398336 0.2827 0.38328
WT Male <-> KO Male 0.2635370 0.2635370 4.2755571 0.1247407 0.0041 0.02130
WT Male <-> KO Female 0.0521060 0.0521060 0.9933835 0.0320515 0.4059 0.40590
WT Female <-> KO Male 0.1359739 0.1359739 2.2360774 0.0693657 0.0607 0.12140
WT Female <-> KO Female 0.0588837 0.0588837 1.1406231 0.0366281 0.3194 0.38328
KO Male <-> KO Female 0.2290033 0.2290033 4.0215307 0.1182055 0.0071 0.02130
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.00865 0.0028848 0.4361    999  0.725
## Residuals 60 0.39688 0.0066146
# Figures
#library packages
pacman::p_load(
  ggthemes,
  ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )

# Genotype

datascores <- as.data.frame(scores(NMDS))# $sites
scores <- cbind(as.data.frame(datascores), Genotype = df_metadata$Genotype)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ Genotype, data = scores, FUN = mean)
seg <- merge(scores, setNames(centroids, c('Genotype','oNMDS1','oNMDS2')),
             by = c('Genotype'), sort = FALSE)

composition_Genotype <- ggplot(scores, aes(x = NMDS1, y = NMDS2, color = Genotype)) + 
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size=0.25) + 
  geom_point(data = centroids, size = 4) +                    
  geom_point()+
  coord_fixed()+ 
  theme_bw() +
  labs(title="Genotype (Distance matrix: Bray-Curtis)")+
  theme(legend.position = "right",
        legend.title= element_text(size=10,face="plain"),
        legend.text = element_text(size=10),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
  scale_y_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
  scale_color_manual(values=c("#3C5488B2","#7E6148B2"))+
  scale_fill_manual(values=c("#3C5488B2","#7E6148B2"))+
  annotate("text",x=-0.5,y=0.75, #0.4 if legend position is right
           label="Stress = 0.14",hjust=0, size=3)+ #2.2 if legend position is right
  annotate("text",x=-0.5,y=0.69,
           label="Pseudo-F = 2.74",hjust=0, size=3)+
  annotate("text",x=-0.5,y=0.63,
           label="p = 0.028",hjust=0, size=3)

# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
        separate(GenotypeSex,c("Genotype","Sex"),sep="_",remove=FALSE)


centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
      separate(GenotypeSex,c("Genotype","Sex"),sep="_",remove=FALSE)

seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
             by = c('GenotypeSex','Genotype','Sex'), sort = FALSE) 


composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size = 0.25) +
  geom_point(data = centroids, size = 4) +
  geom_point()+ #geo,_point()
  coord_fixed()+
  theme_bw() +
  labs(title="")+
  theme(legend.position = "bottom",
        legend.title= element_text(size=10, face="plain"),
        legend.text = element_text(size=14),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 9), size=16),
        axis.title.y=element_text(margin = margin(r = 3), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
  scale_y_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
  scale_color_manual(values=c("#30436d","#9ea5c2","#654e3a","#bcaba6"))+
  scale_fill_manual(values=c("#30436d","#9ea5c2","#654e3a","#bcaba6"))+
  scale_shape_manual(values=c(24,21,24,21))+
  annotate("text",x=-0.5,y=-0.5,
           label="Stress = 0.14",hjust=0, size=4)+
  annotate("text",x=-0.5,y=0.75,
           label=expression(paste("Male (WT vs KO): ",italic("P") == 0.021)),hjust=0, size=4)+
  annotate("text",x=-0.5,y=0.65,
           label=expression(paste("Female (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)


# Print
composition_GenotypeSex_final

ANOISM: Variance between groups and within groups

  • ANOSIM analysis to test whether the variance between groups is larger than within groups (if it is larger, then P values is less than 0.05)
  • Anoism analysis (R value from -1 to 1, R > 0 indicates between-group difference > within group difference; R < 0 indicates within group > between group difference; P value indicates the statistical significance)
# Bray-curtis
dist_bray <- vegdist(df_abundance,  method = "bray")
## Genotype
anosim.genotype_bray <- with(df_metadata, anosim(dist_bray, Genotype))
summary(anosim.genotype_bray)
## 
## Call:
## anosim(x = dist_bray, grouping = Genotype) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.02295 
##       Significance: 0.111 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0260 0.0407 0.0556 0.0681 
## 
## Dissimilarity ranks between and within classes:
##         0%    25%    50%     75% 100%    N
## Between  2 515.75 1035.0 1513.50 2016 1024
## WT       3 506.75  944.0 1453.25 2015  496
## KO       1 481.00 1028.5 1596.25 2011  496
plot(anosim.genotype_bray)

## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)
## 
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.0745 
##       Significance: 0.014 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0333 0.0449 0.0570 0.0770 
## 
## Dissimilarity ranks between and within classes:
##           0%    25%    50%     75% 100%    N
## Between    2 532.75 1031.5 1521.25 2016 1536
## WT Male    3 349.00  784.5 1512.25 2010  120
## WT Female  7 619.75  969.5 1431.00 1934  120
## KO Male    5 457.50 1154.0 1704.25 1979  120
## KO Female  1 395.00  816.5 1192.25 1965  120
plot(anosim.GenotypeSex_bray)